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StellaR: a Package to Manage Stellar 
Evolution Tracks and Isochrones 



by Matteo DeU'Omodarme, and Giada Valle 

Abstract We present the R package stellaR, 
which is designed to access and manipulate pub- 
licly available stellar evolutionary tracks and 
isochrones from the Pisa low-mass database. 
The procedures of the extraction of impor- 
tant stages in the evolution of a star from the 
database, of the isochrones construction from 
stellar tracks and of the interpolation among 
tracks are discussed and demonstrated. 

Due to the advance in the instrumentation, nowa- 
days astronomers can deal with a huge amount of 
high-quality observational data. In the last decade 
impressive improvements of spectroscopic and pho- 
tometric observational capabilities made available 
data which stimulated the research in the globu- 
lar clusters field. The theoretical effort of recov- 
ering the evolutionary history of the cluster bene- 
fits by the computation of extensive databases of 
stellar tracks and isochrones, such as Pietrinferni 
et al. (2006); Dotter et al. (2008); Bertelli et al. (2008). 
We recently computed a large data set of stellar 
tracks and isochrones, "The Pisa low-mass database" 
(DeU'Omodarme et al., 2012), with up to date phys- 
ical and chemical inputs, and made available all 
the calculations to the astrophysical community at 
the Centre de Donnees astronomiques de Strasbourg 
(CDS)^, a data center dedicated to the collection and 
worldwide distribution of astronomical data. 

In most databases, the management of the in- 
formation and the extraction of the relevant evolu- 
tionary properties from libraries of tracks and /or 
isochrones is the responsibility of the end users. 
Due to its extensive capabilities of data manipula- 
tion and analysis, however, R is an ideal choice for 
these tasks. Nevertheless R is not yet well known 
in astrophysics; up to December 2012 only seven as- 
tronomical or astrophysical-oriented packages have 
been published on CRAN (see the taskview Chemo- 
metrics and Computational Physics). 

The package stellaR (DeU'Omodarme and Valle, 
2012) is an effort to make available to the astro- 
physical community a basic tool-set with the follow- 
ing capabilities: retrieve the required calculations 
from CDS; plot the information in a suitable form; 
construct by interpolation tracks or isochrones of 
compositions different to the ones available in the 
database; construct isochrones for age not included 
in the database; extract relevant evolutionary points 
from tracks or isochrones. 

Via anonymous ftp to cdsarc.u-strasbg.fr, or via http: //cdsarc. 



Get stellar evolutionary data 

The Pisa low-mass database contains computations 
classified according to four parameters: the metal- 
licity z of the star, its initial helium value y, the 
value of a-enhancement of the heavy elements mix- 
ture with respect to the reference mixture and the 
mixing length parameter tx.-^i used to model external 
convection efficiency. The values of the parameters 
available in the database can be displayed using the 
function showComposition ( ) : 

> showComposition ( ) 
Mixing-length values: 
1.7, 1.8, 1.9 

alpha-enhancement values: 

0, 1 (i.e. [alpha/Fe] = 0.0 
[alpha/Fe] = 0.3) 



Chemical compositions: 
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The table of chemical composition presents all the y 
values available for a given z. For a set of parame- 
ters, the tracks files are identified specifying the mass 
of the desired model (in the range [0.30, 1.10] Mq 
{Mq = 1.99- 10^3 g is the mass of the Sun), step of 
0.05 Mq), while the age (in the range [8.0 - 15.0] Gyr, 
step of 0.5 Gyr) is required for the isochrones. 

Upon specification of the aforementioned param- 
eters, the stellaR package can import data from CDS 
(via anonymous ftp) from an active Internet con- 
nection. The CDS data are stored in ASCII format, 
and include a header with calculation metadata, such 

-strasbg.fr/viz-bin/qcat?J/A+A/54 0/A26 
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as the metallicity, the inital helium abundance, the 
mixing-length. The import is done via a read. table 
call, skipping the header of the files. 

The following data objects can be downloaded 
from the database site: 

• Stellar track: a stellar evolutionary track com- 
puted starting from Pre-Main Sequence (PMS) 
and ending at the onset of helium flash (for 
masses M > 0.55 Mq) or at the exhaustion of 
central hydrogen (for 0.30 M© < M < 0.50 Mq). 
The functions getTrkO and getXrkSetO can 
be used to access such data; they respectively 
return objects of classes "trk" and "trkset". 

• Stellar ZAHB: Zero-Age Horizontal-Branch 
models. The functions getZahb ( ) can be used 
to access such data; it returns an object of class 
"zahb". 

• HB models: computed from ZAHB to the onset 
of thermal pulses. The functions getHbO and 
getHbgrid ( ) can be used to access such data; 
they respectively return objects of classes "hb" 
and "hbgrid". 

• Stellar isochrones: computed in the age range 
[8 - 15] Gyr. The functions getlsoO and 
getlsoSet () can be used to access such data; 
they respectively return objects of classes "iso" 
and "isoset". 

Readers interested in details about the computa- 
tion procedure are referred to Dell'Omodarme et al. 
(2012). The data gathered from CDS are organized 
into objects of appropriate classes. The package in- 
cludes print and plot S3 methods for the classes 
"trk", "trkset", "hb", "zahb", "hbgrid", "iso", and 
"isoset". 

As an example, we illustrate the recovering of the 
stellar track for a model of mass M = 0.80 Mq, metal- 
licity z = 0.001, initial helium abundance y = 0.25, 
mixing-length = 1.90, a-enhancement [a/Fe] = 
0.0. 



> track <- getTrk(m =0.80, z 

y = 0.25, ml = 1.9C 

> track 
Stellar track 

Mass =0.8 Msun 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 

> names (track) 

[1] "mass" "z" "i 

"alpha. enh" "data" 

> class (track) 

[1] "trk" "stellar" 



0.001, 

afe = 0) 



"ml" 



The function getTrkO returns an object of class 
"trk", which is a list containing the track meta- 
data, i.e. as the star mass, the metallicity, the ini- 
tial He abundance, the mixing-length and the a- 
enhancement, and the computed data in the data 
frame data. Track data contains the values of 15 vari- 
ables: 



> names (track$data) 
[1] "mod" "time" 
"mass" "He" 
"MHEc" "Lpp" 



"logL" "logTe" 
"logic" "logRHOc" 
"LCNO" "L3a" 



"Lg" 



"radius" "logg" 



The included variables are: mod the progressive 
model number; time the logarithm of the stellar age 
(in yr); logL the logarithm of the surface luminos- 
ity (in unit of solar luminosity); logTe the logarithm 
of the effective temperature (in K); mass the stellar 
mass (in unit of solar mass); He the central hydrogen 
abundance (after hydrogen exhaustion: central he- 
lium abundance); logic the logarithm of the central 
temperature (in K); logRHOc the logarithm of the cen- 
tral density (in g/ cm^); MHEc the mass of the helium 
core (in units of solar mass); Lpp the luminosity of 
pp chain (in units of surface luminosity); LCNO the lu- 
minosity of CNO chain (in units of surface luminos- 
ity); L3a the luminosity of triple-a burning (in units 
of surface luminosity); Lg luminosity of the gravita- 
tional energy (in units of surface luminosity); radius 
the stellar radius (in units of solar radius); logg the 
logarithm of surface gravity (in cm/ s^). 

Similarly the part of the track starting from ZAHB 
and ending at the onset of thermal pulses can be 
downloaded by the call: 



> hbtk <- getHb(m = 0.80, 

y = 0.25, ml 

> hbtk 

Stellar track from ZAHB 

Mass =0.8 Msun 
Mass RGB =0.8 Msun 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 



z = 0.001, 
1.90, afe = 0) 



> names (hbtk) 
[1] "mass" 
"ml" 



"massRGB" "z" 
"alpha. enh" "data" 



> class (hbtk) 

[1] "hb" "stellar" 

which returns an object of class "hb", which differs 
from an object of class "trk" only for the presence of 
the variable mas sRGB, i.e. the Red-Giant Branch (RGB) 
progenitor mass. 

Usually a set of tracks with different mass 
and/ or metallicity are needed for computations. The 
package stellaR provides the function getlrkSetO, 
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which can download a set of tracks with differ- 
ent values of mass, metallicity, helium and mixing- 
length. As an example the whole set of masses (from 
0.30 to 1.10 Mq, step of 0.05 M©), for metallicity 
z = 0.001, initial helium abundance y = 0.25, mixing- 
length a^i = 1.90, a-enhancement [a/Fe] = 0.0 can be 
downloaded as follows: 

> mass <- seq(0.3, 1.1, by = 0.05) 

> trks <- getTrkSet (m = mass, z = 0.001, 

y = 0.25, ml = 1.90, afe = 0) 

> trks 
[[1]] 

Stellar track 

Mass =0.3 Msun 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 



[[2]] 

Stellar track 

Mass =0.35 Msun 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 




3.90 3.85 3.80 



3.75 3.70 
log Teff 



3.65 3.60 3.55 



Figure 1: The evolutionary tracks for masses from 
M = 0.30 M© to M = 1.10 Mq from PMS to He flash. 
The parameters of the calculations are: z = 0.001, y = 
0.25, ajni = 1.90, [a/Fe] =0.0. 

The function getTrkSet () returns an object of 
class "trkset", a list containing objects of class 
"trk". The track set can be displayed in the usual 
(logTgff, log L/L©) plane by a call of the function 

plotO: 



> plot (trks, Ity = 1:2) 

The output of the function is shown in Figure 1. The 
plot is produced by a call to the function plotAstro ( ) 
which allows the user to customize several aspects of 
the plot, such as the axes label, the number of minor 
ticks between two major ticks, the limits of the axes, 
the color and type of the lines (as in the example), the 
type of the plot (lines, points, both, . . . ). 

The use of plotO function is further demon- 
strated in Figure 2, where, for z = 0.001, y = 0.25, a^i 
= 1.90, [a/Fe] = 0.0, are displayed the evolutionary 
track for M = 0.80 M© from PMS to He flash (black 
line) and from ZAHB to thermal pulses (green line). 
The figure is obtained as follows: 

> plot (track) 

> plot(hbtk, add = TRUE, col = "green") 




3.85 3.80 3.75 3.70 3.65 
log Teff 



3.60 



3.55 



Figure 2: (black line): evolutionary tracks for mass 
M = 0.80 Mq, z = 0.001, y = 0.25, a^i = 1.90, [a/Fe] = 
0.0. from PMS to He flash, (green line): evolutionary 
track from ZAHB to thermal pulses. 

Apart from the plots discussed before, it is easy 
to display other relation between the computed vari- 
ables. In the following example we get the data for 
two masses, namely 0.50 and 1.00 M© and plot the 
trend of the radius (in units of solar radius) versus 
the logarithm of the age for the first 100 models. The 
resulting plot is displayed in Fig. 3. 

> trkr <- getTrkSet (m = c(0.5, 1), z = 0.01, 

y = 0.25, ml = 1.8, afe = 0) 

> mydata <- do. call (rbind, lapply(trkr, 

"[[", "data")) 

> D <- mydata [mydata$mod <= 100, ] 

> key <- as . numeric (factor (D$mass) ) 

> plotAstro (D$time, D$radius, type = "p", 

pch = key, ylab = "Radius (Rsun) ", 
xlab = "log age (yr)") 
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> legendC'topright", c("M=0.50", "M=1.00"), 
pch = 1:2) 
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4 5 6 7 8 9 

log age (yr) 

Figure 3: Radius versus the logarithm of the age for 
the first 100 models of two stars with different mass 
(M = 0.50 Mq, and M = 1.00 M©) and identical com- 
position, z = 0.01, y = 0.25, ftmi = 1.8, [a/Fe] = 0.0. 

Isochrones can be obtained from CDS and plotted 
in a similar way. As an example, let get isochrones of 
9 and 12 Gyr for z = 0.001, y = 0.25, a^l = 1-90, [a/Fe] 
= 0.0: 

> isc <- getlsoSet (age = c(9, 12), z = 0.001, 

y = 0.25, ml = 1.90, afe = 0) 

> isc 
[[1]] 

Stellar isochrone 

Age = 9 Gyr 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 

[[2]] 

Stellar isochrone 

Age =12 Gyr 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 

attr(, "class") 

[1] "isoset" "stellar" 



> names (isc [ [1] ] ) 
[1] "age" "z" 

"alpha. enh" "data" 

> names (isc [ [1] ] $data) 



"ml" 



[1] "logL" "logTe" "mass" "radius" "logg" 

The function returns an object of class "isoset", a list 
contarnrng objects of class "iso". These last objects 
are lists containing metadata (age, metallicity initial 
He abundance, mixing-length, and a-enhanchment) 
and the data frame data, which contains the com- 
puted theoretical isochrones data. Figure 4 shows the 
set of isochrones plotted with the commands: 

> plot (isc, Ity = 1:2) 

> legend("topleft", c("9 Gyr", "12 Gyr"), 

Ity = 1:2) 



^ -| 9 Gyr 

12 Gyr 
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Figure 4: Isochrones in the theoretical plane for 9 
and 12 Gyr. Calculations performed for z = 0.001, y = 
0.25, a^l = 1-90, [a/Fe] = 0.0. 

Tools for interpolating among data 
structures 

Even if the database provides a fine grid of mod- 
els it is possible that the specific model needed 
for a comparison with observational data has not 
been calculated. To address this issue, the package 
stellaR includes tools for construction of new sets 
of isochrones. The simplest case is whenever one 
desires an isochrone for ages not included in the 
database, but for a combination of metallicity. He 
abundance, mixing-length and a-enhanchment exist- 
ing in the database. The function make I so () is able 
to compute the required isochrones by mean of in- 
terpolation on a tracks set. The interpolation can be 
performed for age in the range [7 - 15] Gyr. The user 
has the choice to explicitly provide to the function a 
set of tracks, previously downloaded from CDS, or 
to specify the required composition of the tracks to 
be downloaded for the interpolation. To show the 
usage of the function we use the object trks down- 
loaded before to obtain an isochrone of age 9.7 Gyr: 
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> iso.ip <- makeIso(age = 9.7, tr = trks) 

> iso.ip 

Stellar isochrone 

Age = 9.7 Gyr 
Z = 0.001 , Y = 0.25 
Mixing length = 1.9 
[alpha/Fe] = 

the call produces a result identical to make I so (age 
= 9.7, z = 0.001,y = 0.25, ml = 1.9,afe = 0); in 
this last case the data are taken from CDS before the 
interpolation procedure. 

The interpolation technique is based upon the 
fact that all the tracks contain the same number of 
points by construction, and that a given point cor- 
responds to the same evolutionary phase on all the 
tracks. We define S(M) as the set of tracks to be used 
for interpolation, parametrized by the value of the 
mass M. Let f,(M) be the evolutionary time for the 
fth point on the track of mass M, and A be the age of 
the required isochrone. Let k be the point on the track 
of lower mass of the set S(M) for which tk{M) > A. 
For each point > A: on the tracks in S(M), a linear 
interpolation in age of the values of mass, logarithm 
of the effective temperature and logarithm of the lu- 
minosity is performed among tracks. These points 
define the required isochrone. A potential problem 
of this simple procedure will occur whenever mas- 
sive stars develop a convective core during the Main 
Sequence (MS). In this case, as shown for example in 
Mowlavi et al. (2012), the monotonic trend of the evo- 
lutionary time - that decreases with increasing stellar 
mass at the end of the MS - inverts at the middle of 
the MS. However the problem will be encountered 
only for early -time isochrones, for which the mass at 
the isochrone Turn-Off will be in the interval during 
which the convective core develops. The procedure 
outlined in this Section is adequate for construction 
of isochrones throughout the range allowed by the 
function. 

Tracks interpolation 

The package stellaR provides also a tool for perform- 
ing a 3D interpolation on the database to construct a 
set of tracks for values of metallicity, helium abun- 
dance and mixing-length not included in the compu- 
tations available at CDS. The function interpTrkO 
can be used for this procedure. A call to this function 
causes the download from CDS of the sets of tracks 
needed for the interpolation. 

The new set of tracks is computed by mean 
of a linear interpolation. The metallicity is log- 
transformed before the interpolation procedure. Let 
Tz y (M), be fth point on the dataset containing the 
evolutionary time, the effective temperature and the 
logarithm of surface luminosity for the track of mass 
M, and given composition. The interpolation algo- 



rithm proceeds as follows: 

8 sets 4 sets 2 sets 1 sets 

where the symbol * means that interpolation oc- 
curred in the substituted variable. The selection 
of the tracks set which enter in the interpolation is 
based upon the identification of the vertexes of the 
cell of the (z, y, a^i) space containing the point iden- 
tified by the required parameters. Then, for all the 
17 masses at the vertexes, the linear interpolation 
described above is performed. In the worst case 
scenario, whenever no one of the supplied parame- 
ters values exist in the database, the interpolation re- 
quires 23 =8 sets of 17 masses. The algorithm is how- 
ever able to reduce the dimensionality of the process 
if some of the variables values exist in the database. 

As a demonstration, let us compute a set of tracks 
with mixing-length value = 1.74, z = 0.002, y = 
0.25, [a/Fe] = 0.0: 

> ip.trk <- interpTrk(z = 0.002, y = 0.25, 

ml = 1.74, afe = 0) 

Since the values of z and y exist in the database, only 
an interpolation on the mixing-length value is per- 
formed by the code. The set of tracks can be used for 
isochrones construction, like a standard set of tracks: 

> ip.iso <- makeIso(age = 12, tr = ip.trk) 

Keypoints extraction 

Important stages in the evolution of a star are defined 
as "keypoints", e.g. hydrogen core exhuastion, Turn- 
Off luminosity, RGB tip luminosity. To simplify their 
extraction the package stellaR provides the function 
keypoints 0, which operates on an object of class 
"trk" or "iso". 

The function extracts from the data stored in ob- 
jects of class "trk" the rows of the data frame relative 
to the following evolutionary stages: 

1. ZAMS. Zero-Age Main-Sequence, defined as 
the point for which the central H abimdance 
drops below 99% of its initial value. 

2. TO. Turn-Off, defined as the point for which 
the effective temperature reaches its maximum 
value. If multiple lines satisfy the constraint, 
the values of all the rows are averaged. 

3. BTO. Brighter Turn-Off, defined as the point for 
which the effective temperature drops below 
the one of the TO minus 100 K. The points can 
not exist for low masses. Details on the advan- 
tages of this evolutionary point with respect to 
the TO can be foimd in Chaboyer et al. (1996). 
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4. exHc. Central H exhaustion, defined as the 
point for which the central H abundance is 
zero. For low masses the point can coincide 
with TO. This is the last point of the tracks with 
mass lower or equal to 0.50 Mq. 

5. Heflash. Helium flash, the last point of the 
track for masses higher than 0.50 Mq. 

When the function is called on an object of class 
" i s " it returns a data frame containing only TO and 
BTO phases. 

In both cases the function inserts in the re- 
turned data frame the columns relative to mass 
(or age for isochrones), metallicity, initial He value, 
mixing-length, a-enhancement, and evolutionary 
phase identifier. 

As a demonstration we extract the TO and BTO 
points from the object isoc generated in a previous 
example: 



> kp <- keypoints (isoc) 

> kp 

logL logTe 
0.4044723 3.816652 
0.5556971 3.809943 
0.2917250 3.803611 
BTOl 0.4495597 3.796545 
logg age z 
4.178911 
4.009265 
4.205965 
BTOl 4.027426 
id 

TO 1 
BTO 2 
TOl 1 
BTOl 2 



The points can be easily superimposed to the 
isochrones in the theoretical plane. The top panel of 
Figure 5 is obtained with the following commands: 



TO 

BTO 

TOl 



TO 

BTO 

TOl 



9 0.001 

9 0.001 

12 0.001 

12 0.001 



mass radius 
,8396903 1.23558 
,8561162 1.51671 
,7772973 1.15234 
,7909500 1.42768 

y ml alpha. enh 
,25 1.9 
,25 1.9 
,25 1.9 
,25 1.9 



> plot (isoc) 

> points (logL 



kp, pch = id, 



logTe, data 
col = id) 

> legendC'topleft", c("TO", "BTO"), 

pch = 1:2, col = 1:2) 

As a last example we extract a set of tracks for masses 
in the range [0.40 - 1.10] M© and three metallicity z 
= 0.0001, 0.001, 0.01 and we display (bottom panel in 
Figure 5) the time of exhaustion of central hydrogen 
as a function of the mass of the star. 

> mT <- seq(0.4, 1.1, by = 0.1) 

> zT <- c(0.0001, 0.001, 0.01) 

> tr <- getTrkSet (m = mT, z = zT, y = 0.25, 

ml = 1.9, afe = 1) 

> kp <- keypoints (tr) 

> kpH <- kp [kp$id == 4, ] 



> xlab <- expression (italic (M) ~ 

(italic (M) [sun] ) ) 

> ylab <- "log age (yr) " 

> symbol <- as .numeric (factor (kpH$z) ) 

> plotAstro (kpH$M, kpH$time, type = "p", 

pch = symbol, xi = 0.1, xlab = xlab, 
ylab = ylab) 
lab <- levels (factor (format (kpH$z, 

nsmall = 4, scientific = FALSE))) 

> legend("topright", lab, pch = 1:3) 



o 



o TO 
A BTO 




3.85 



3.80 



3.75 



3.70 
log Teff 



3.65 



3.60 



3.55 



to 



o 0.0001 ■ 
& 0.0010 
+ 0.0100 ■ 



0.4 



0.5 



0.6 



0.7 0.8 



0.9 



1.0 



1.1 



Figure 5: Top panel: same as Figure 4. The position 
of Turn-Off (TO) and Brighter TO (BTO) are shown. 
Bottom panel: time to exhaustion of central hydro- 
gen as a function of the mass of the star. The symbols 
identify the values of the metallicity z form 0.0001 to 
0.01. 
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Summary 

This paper demonstrated how the package stellaR 
can be useful to the astrophysical community as 
a tool to simplify the access to stellar tracks and 
isochrones calculations available on-line. A set of 
tools to access, manipulate and plot data are in- 
cluded in the package and their usage was shown. 
The interpolation functions included in the package 
can be used to safely produce tracks or isochrones at 
compositions not included in the database, without 
the need that users develop software on their own. A 
planned extension of the package is the modification 
of the algorithm of isochrone construction to make 
feasible the calculation of isochrone of young ages. 
This step can be useful in view of a possible exten- 
sion of the Pisa database to higher masses, or to ma- 
nipulation of data stored in other databases. In fact, 
while the package is currently developed for access- 
ing data from Pisa low-mass database, other public 
databases can be in principle accessed in the same 
way. This step is however complicated by the fact 
that no standard for stellar model output exists in the 
astrophysical community, requiring the personaliza- 
tion of the interface functions for each set of data. 
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